Introduction to spatial models for citizen science data

Before jumping into the extra complexities of spatial models, it is worth reflecting when we do not need spatial models:

1. spatial field not of primary interest to the study AND
2. spatial variation is captured by covariates included in the model

1. Spatial field of interest (distribution models …)
2. Covariates cannot account for the spatial signal

Biodiversity recording from citizen science programms often contain a spatial signal due to citizen reporting biodiversity records in specific areas:

1. Easily accessible (close to roads, close to [vacation] home …)
2. Areas with interesting species / habitats

In addition, citizen science biodiversity monitoring programs give informations on where a species was recorded (and is therefore assumed to be present), but an absence of records cannot be usually assumed to represent a species absence. This type of data is usually called presence-only.

One interesting class of models to fit to presence-only data with spatial signal is: point pattern analysis.

Two main packages to fit point patterns in R:

1. spatstat: not covered in this workshop
2. INLA/inlabru

INLA stands for Integrated Nested Laplace Approximation, it is an approximate Bayesian approach that offers a faster alternative to other Bayesian algorithms such as MCMC or HMC.

INLA can fit a wide range of models including:

1. Hierarchical models (random effect)
2. Autocorrelated models (like AR(1))
3. Spatial models (including point patterns)
4. Space-time models

Spatial effect in INLA are fitted using Stochastic Partial Differential Equations (SPDE)

Picture from Krainski et al.

The continuous spatial effect is estimated by using basis functions (like smooths) on a mesh which dramatically speed up the estimation of continuous spatial effects [more on this later].

Picture from Krainski et al.

Given a mesh, to fit a spatial model in INLA/inlabru we need to define the SPDE and more specifically the priors for the standard deviation (sigma) and the scale (range) of the spatial effect.

There are various ways to define the SPDE, the recommended one is via Penalized Complexity priors that is parametrized as follow:

1. variance of the spatial effect: \( P(\sigma > \sigma_0) = s \), where \( \sigma_0 \) and s are two values to be specified
2. range of the spatial effect: \( P(\rho < \rho_0) = p \), where \( \rho_0 \) and p are two values to be specified

From Fulgstad et al 2019:

In the simulation study we observe good coverage of the equal-tailed 95% credible intervals when the prior satisfies \( P(\sigma > \sigma_0) = 0.05 \) and \( P(\rho < \rho_0) = 0.05 \), where \( \sigma_0 \) is between 2.5 to 40 times the true marginal standard deviation and \( \rho_0 \) is between 0.1 and 0.4 of the true range.

In other words if we assume a standard deviation of the spatial effect to be around 2 and a range of 1 we can set the spde in R:

spde <- inla.spde2.pcmatern(mesh,
        # P(sigma > 10) = 0.05
        prior.sigma = c(10, 0.05),
        # P(range < 0.1) = 0.05
        prior.range = c(0.1, 0.05))

Ok, but what do this sigma and range really represent?

Under the hood INLA uses the Matern correlation / covariance function to estimate the spatial field.

Remember that our spatial models are usually:

\[ y(s_i) \sim D(\mu(s_i), ...) \\ \mu(s_i) = \beta X + u(s_i) \]

Where \( s_i \) are the n locations where the data were sampled and u is a continuous spatial field that is commonly assumed to follow a multivariate normal distribution that can be parametrized as a gaussian process:

\[ u \sim GF(0, \Sigma) \]

The \( \Sigma \) matrix represent the correlation / covariance between the n observations, in other words if you have 100 locations, \( \Sigma \) has 100 rows and 100 columns.

Fortunately we do not have to individually estimate all of these parameters, we can use the Matern function to compute the correlation / covariance between all observations given their (eucledian) distance following:

\[ cov(s_i, s_j) = Matern(\|s_i - s_j\|, \sigma, \rho) \]

Where \( \sigma \) is the variance in the spatial effect and \( \rho \) the range.

mesh <- inla.mesh.2d(loc = data, max.edge = 0.1)
spde <- inla.spde2.pcmatern(mesh,
        prior.sigma = c(10, 0.05),
        prior.range = c(0.1, 0.05))
formula <- coordinates ~ Intercept + 
  spat(coordinates, model = spde)
mod <- lgcp(components = formula,
            data = spdf,
            samplers = spol)

Where spdf is a SpatialPointDataFrame of the observations and spol is a a SpatialPolygonDataFrame of the domain boundary.

range <- spde.posterior(mod, "spat", "range")
sigma <- spde.posterior(mod, "spat", "variance")

matern <- spde.posterior(mod, "spat", "matern.correlation")